Elastic Constants of Quantum Solids by Path Integral Simulations 
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Two methods are proposed to evaluate the second-order elastic constants of quantum mechani- 
cally treated solids. One method is based on path-integral simulations in the NVT ensemble using 
an estimator for elastic constants dj . The other method is based on simulations in the NpT ensem- 
ble exploiting the relationship between strain fluctuations and elastic constants. The strengths and 
weaknesses of the methods are discussed thoroughly. We show how one can reduce statistical and 
systematic errors associated with so-called primitive estimators. The methods are then applied to 
solid argon at atmospheric pressures and solid helium 3 (hep, fee, and bcc) under varying pressures. 
Good agreement with available experimental data on elastic constants is found for ''He. Predictions 
are made for the thermal expectation value of the kinetic energy of solid ^He. 
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I. INTRODUCTION 



The quantum nature of atomic motion alters thermal 
and mechanical properties of condensed matter at low 
temperatures. The most striking feature is the well- 
known quantum mechanical freezing of solids below the 
Debye temperature: Specific heat and thermal expansion 
vanish as the temperature approaches absolute zeroEJ. 
This behavior is different from what we would expect 
from classical statistical mechanics. A "classical solid" 
typically has a positive expansion coefficient at low tem- 
peratures and a specific heat close to Sfes per atom. 

Path integral Monte Carlo (PIMC)ErEl and path inte- 
gral molecular dynamics (PIMD)LfLl are convenient nu- 
merical methods to predict the proper low-temperature 
properties of condensed matter, provided that model po- 
tentials are available that describe the interatomic inter- 
actions sufficiently well. So far, most PIMC and PIMD 
studies have been restricted to the calculation of struc- 
tural and thermal properties of quantum solids or to the 
calculation of equations of state of condensed rare gases. 
The computation of the tensor of the elastic constants, 
which is an important material property in many engi- 
neering applications, however, has not been put forward 
yet in the framework of path-integral simulations. 

The proper evaluation of the elastic constants is far 
more complicated than textbooks on solid state physics 
make us believe, because it is not sufficient to calculate 
second derivates of interaction potentials at equilibrium 
(or thermal) positions and construct from thiS|_the vari- 
ous Cij. As shown by Squire, Holt, and Hoover J3 thermal 
fluctuations of a generalized strain tensor can alter elastic 
constants significantly. Of course, quantum fluctuations 
can be expected to result in a similar effect. These fluctu- 
ations vanish at small temperatures for classical systems, 
but they do not necessarily vanish for quantum mechani- 
cal systems. Another difficulty is that quantum solids do 
not adopt a configuration at T = K where the poten- 
tial energy surface is minimized. Instead, the atoms also 



probe non-harmonic parts of the potentials due to quan- 
tum fluctuations, typically leading to lattice constants 
that are slightly increased with respect to the equivalent 
classical system. 

In this paper, we want to propose two methods which 
enable us to determine accurately elastic constants of 
solids such that quantum effects of the atomic motion 
are fully included in the treatment. One possibility is 
to use the definition of the elastic constants as the sec- 
ond derivative of the free energy with respect to strain 
tensor elements and evaluate the final expression in the 
NVT ensemble. This procedure is similar to the method 
proposed by Squire, Holt, and Hoover,Q with the differ- 
ence that our partition function is quantum mechanical. 
Alternatively, one can perform simulations in the NpT 
ensemble and relate the fluctuations of the strain tensor 
to the elastic constants as has been originally done by 
Parrinello and Rahman for classical systems Jj 

Both methods are prone to produce large statistical 
error bars. In classical simulations, the computation of 
elastic constants using the NpT ensembles is knpvn to be 
much less efficient than in the NVT ensemblcEj. How- 
ever, the estimator derived for the quantum simulations 
in the NVT ensemble can be expected to become un- 
reliable when the Trotter number P in the path-integral 
simulation becomes large, just like the so-called primitive 
estimator for the kinetic energyO. It is thus necessary to 
discuss the prospective statistical errors in detail. 

It is important to note that presently only isothermal 
quantum mechanical elastic constants can be calculated. 
Adiabatic elastic constants can only be obtained if simu- 
lations are.carried out in the NVE ensemble or the NpH 
ensemble.EI Constraining the energy E or the enthalpy H 
to constant values in the simulation does not necessar- 
ily correspond to conserved E or H oi the real quantum 
system. 

The remainder of this paper is organized as follows: In 
Sec. H, the equations are derived that allow the correct 
determination of all Cij in the NVT ensemble in terms of 
a path integral formulation. In this context, we will give 



an in depth discussion of the statistical errors that will 
arise as a consequence of the corrections to Cy associ- 
ated with the quantum mechanical kinetic energy. It will 
be shown, how, simple corrections to primitive estimators 
can reduce their statistical incertainties considerably. We 
will also briefly review the method to determine Cij in 
the NpT ensemble as well as the simulation algorithm 
that is used for the applications. In Sec. Ill, the meth- 
ods will be applied to the calculation of elastic constants 
of solid Argon as well as of fee, bcc, and, hep '^He. In the 
case of '^He, the corrections to primitive estimators will 
also be used to make predictions for thermal expectation 
values of the kinetic energy. Conclusions are drawn in 
Sec. fV, 



II. THEORY 

A. Derivation of Elastic Constants 

As pointed out corpctly for classical systemsB, isother- 
mal elastic constantslla Cy are insufficiently described if 
they are evaluated solely on the basis of averageing the 
so-called Born contribution C, 
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which is the thermal average of the second derivative of 
the potential energy Vpot with respect to strain tensor 
elements e^ and ej in the NVT ensemble. The full elastic 
constants are obtained if (V^ot) is replaced with the free 
energy F{N,V,T) = -kBTlnZ{N,V,T), e.g.. 
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with Z{N, V, T) the isothermal partition function. This 
proper definition results in corrections to C 
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leading correction terms C^ ■ are fluctuations of the 
(generalized) instantaneous strains, 
the form 
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Note that one has to properly symmetrize when using the 



Voigt notation, e.g., C 
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contribution can usually be neglected in classical simu- 
lations, knowing the correct form is crucial for quantum 
systems as can be seen further below. Thus, in classical 
solids at non-zero temperaures, the elastic constants can 
be estimated as C,, = c{f ""^ + C^^""'^ + cj'"\ 

The quantum statistical formulation of the partition 
function of A^ identical particles with mass m in terms of 
a path-integral formulationtj is given by 
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with A(/3/P) the thermal de Broglie wavelength at tem- 
perature f3/P and the effective potential Ves 



N P 



^- = EE 



ra=l t=l 



rn I Rit — Rjt+i 

' [ i3h/p 



7 ^ V'pot \Rnn' t) 



(7) 



Here, the coordinate of particle n at "Trotter time" t is 
denoted as Rnt and Rnn' t denotes the distance between 
particle n and n' at "Trotter time" t. If exchange effects 
are neglected, one quantum point particle is represented 
as a closed classical ring polymer with the boundary con- 
dition Rnt = Rnt+p- The temperature at which the sim- 
ulation of the representation of the quantum particles is 
done is PT. Strictly speaking, the quantum limit is only 
obtained for infinite large Trotter numbers P, however, 
for practical purposes, it is usually sufficient to choose 
PT in the order of two times the Debye temperature. 

In the following, it is convenient to represent the po- 
sition of the particles in reduced dimensionless variables 
fnt and a (symmetric) matrix hafs that contains the shape 
and the size of the simulation cell such that 



where again the expectation values are evaluated in the 
NVT ensemble. As far as classical elastic constants are 
concerned, the only missing contributions C^^ '" to the 
correct Cij stem from the ideal gas part of the partition 
function. These kinetic corrections are given by 
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or in tensor notation for cubic symmetry 
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with < rnta < 1 and V — det h. Spatial periodic 
boundary conditions are applied by subtracting or adding 
unity to Vnta once a molecular dynamics step has moved 
fnta out of the allowed range. With this transformation, 
the integration over J d'^Rn ■ ■ ■ J d'^R^p in Eq. (0) can 
be replaced with the expression V^^ J d'^ru ■ ■ ■ J drrj^p. 



This makes it possible to take the derivative of 
In Z{N,V,T) with respect to the strain e^. If we do not 
use the Voigt notation, a (virtual) variation in the stress 
tensor Ssa/s can be expressed as 
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Combining Eq. (||) with Eqs. (||) through (^) then results 
in contributions to the elastic constants that resemble 
those obtained for classical systems. In particular, l^pot 
in Eqs. @j and (y|) has to be replaced with Vcs/ P, and 
the factor NksT in Eqs. |j and || has to be replaced with 
NksTP. 

If elastic constants are evaluated in an NpT ensemble 
use can be made of the relation 



{5e^p5ep,) - {ksT/V) (C-i)^^_^, 
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where C has to be represented as a 6 x 6 matrix, where 
we again return to the Voigt notation. 



B. Application to Ideal Gas 

In order to discuss the expected statistical errors, it is 
helpful to consider different contributions to the elastic 
constants. In order to do this, we split the effective po- 
tential T4ff into two parts, the real potential l^ot and 
the remaining part of the right hand side of Eq. (1^), 
which we call Vq. The terms contributing to Cij can 

then be decomposed into the Born contribution C\^ ° , 

a term Q^ om-qj^ which is obtained by replacing l^ot 

in the definition of C,;, °'" with Vq, the term associated 

with the fluctuations of the generalized stress C} "'^ , a 

similarly defined term C^- "'^ stemming from fluctua- 



tions of dVq/dci, the cross correlation C, 
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All terms related to Vq increase linearly in P in leading 
order as one can show in a particularly compact form if 
we exploit the harmonic character of the springs in the 
"ring polymers" , by introducing appropriate normal co- 
ordinates 
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so that V^ is diagonal in B 
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with fcg = 4msin2(7rg/F) / {PTi/Pf 
For the ideal gas we obtain: 
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while the estimator in the presence of a non-vanishing 
potential have the form 
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(18) 



^^ST'^ = -2NkBT{P ~ l)S^^dps, (19) 

with oBnta — Rnta — Rnt+la- 

In the presence of a non- vanishing V^ot, the individ- 
ual terms will only be effected slightly, because in the 
quantum limit, changes in the interaction potential are 
much smoother than the energy changes due to changes 
in Vq. The net and properly symmetrized Cij, however, 
do have well-defined values in the limit P -^ oo, which 
are, of course, sensitive to the interaction potential Vpot • 
For the ideal gas, the symmetrized results are listed in 
Table fl In a numerical calculation, it clearly will be a 
problem that individual terms in Eqs. ( [l4| - |lq ) are of order 
P > > 1 while the final result is of order unity. 



C. Improved Primitive Estimators in PIMD 

Here, we want to suggest how one can easily improve 
on the statistical properties of so-called primitive esti- 
mators in PIMD simulations. As a first step we define 
a function 6Kap which measures how far the averaged 
tensor of the kinetic energy deviates from its expectation 
value 
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n,t 



■nqct 



where rhit represents the (kinetic) mass of particle n at 
Trotter time t (note that kinetic masses carL-jbe chosen 
arbitrarily independent of the physical masaj) and Via 
is the a component of its velocity. (•)simui. symbolizes 
an expectation value obtained in a simulation. Ideally, 
5Kap would be zero, but, due to finite statistics, finite 
time-steps, and other round-off errors, we will usually 
find SKap ^ 0. It is very likely that this deviation SKap 
will convert into a similarly large deviation in potential 
energy, in particular into Vq for large Trotter numbers P. 



At large P, the external potential is locally only a small 
perturbation to the springs connecting neigbored beads 
on the ring. Due to the conversion of kinetic energy to 
potential energy, similar deviations from the exact ther- 
modynamic average can be expected in Vq. This makes 
it possible to define an optimized estimator for quantities 
associated with V^, namely 
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Similar correction terms can easily be generalized to 
higher orders, e.g. deviatians of {rnffViaViisVi^Vis) simu. 
from its thermal expectation value can be expected to 
convert into {{mP'^ / (3h)^SRiaSRii3SRi^SRis)simu. 

In passing we want to explicitly give the improved esti- 
mator for the kinetic energy, which can be usedin PIMD. 
It is similar to the so-called primitive estimatoilij, but the 
expression iNksTP/i arising in the original estimator is 
replaced by the actual average net kinetic energy. With 
the quantities introduced in Eq. rtO), the average kinetic 
energy may simply be expressed as 
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As outlined above, the new primitive estimator benefits 
from the correlation of kinetic and potential energy in 
the isomorphic classical representation. This is demon- 
strated in Fig. n^, where the average statistical deviation 
from terms of the type Cim are shown. 




100 



Note that the statistical error of AC}]^"'^" , which is 
identical to the statistical error of the original primitve 
estimator, increases much slower than the statistical pt.- 
ror of the primitve estimator when examined originallyllil. 
This is due to the dcYf lopment of efficient sampling meth- 
ods in the meantimeEl which suppress slowing down with 
increasing P. The new estimator has a strongly reduced 
statistical noise as compared to the original one, yet, re- 
sults in the correct expectation value. Furhtermore, using 
larger time steps than in our production runs, we have 
noticed that systematic errors due to finite-time steps, 
are considerably reduced as well with the new estima- 
tor. We want to note that in later production runs for 
strongly quantum mechanical systems such as solid ^He, 
statistical error bars obtained with the new estimators 
are as small as statisticai-error bars asscociated with the 
so-called virial estimatoiO. 



D. Elastic Constants Under Pressure 

Care has to be taken when discussing elastic constants 
of systems under pressure, because their definition is not 
unique.Ej One commonly distinguishes between the so- 
called Birch coefficients Bij and the elastic constants Cy . 
They are defined as the second derivative of the free en- 
ergy with respect to the Lagrangian strain and the "reg- 
ular" strain, respectively. The relation between Bij and 
Cij merely depend on the externally applied stress but 
not on the symmetry of the crystal. In thei-Gase of hy- 
drostatic pressure, the relations are given bytJ 



Bi 
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Ay = 



(23) 
-p for z < 3 and 



with Aij = —p for 1 < i < 6, 
j < 3 with i ^ j , and zero else. 

The generalization of our previos estimator of Cij to 
the non-zero pressure case are as follows: Evaluation of 
the strain fluctuations as indicated in Eq. 11^ result in 
the Birch coefficients, while the formula given for zero- 
pressure Cij in the NVT ensemble enable us to calculate 
the Cij at non-zero pressures. 

Note that the proper stability criterion for solids un- 
der pressure is a positive definite matrix of Birch coeffi- 
cients tather than a positive definite matrix of elastic con- 
stants.Ej It is important to keep in mind that symmetry 
relations of Cij which are obtained under the assumption 
of short-ranged two-body potentials are not valid for^B^ 
in the case of any non-zero externally applied stress.Ej 



E. Simulation Method 



FIG. 1. Statistical variances of the components C''^'"' and 
C °''"~'i' for varios Trotter numbers P. Similar conditions 
like same number of MD steps, same temperature, etc. have 
been used for all calculations. 



Simulations are carried out in both, the NpT and the 
NVT ensemble. Molecular dynamics are favorable over 
Monte Carlo simulations in the constant stress ensem- 
ble {NpT), because all variables are averaged simulta- 
neously - in particular all elements of the strain tensor. 



In order to keep the externally applied-stress tensor con- 
stant, the Parrinello-Rahman methodllj has been applied 
to the classical representation which is isomorphic to the 
quantum system. The classical representation is defined, 
by Eqs. (ph and (R). More details are given elsewheretZl 
on how to efficiently collapse the time-scales associated 
with the different motions of the system such as the intra- 
molecular breathing of the closed classical ring represent- 
ing the quantum point particle, the center-of-mass mo- 
tion of the ring, and the flucuations of the cell geometry. 



III. APPLICATIONS 
A. Argon 

Argon is a convenient test case for the calculation of 
quantum-mechanical elastic constants, because it lies in 
between what is considered a quantum solid like solid he- 
lium under pressure, and what is considered a "classical" 
solid like solid Xenon. For the Argon test-case, we intend 
to determine the shift of the quantum mechanical elastic 
constants with respect to the classical elastic constants 
and the relative importance of the individual contribu- 
tions to the net result. The results will be an indicator 
for what we can expect from the quantum mechanical 
shift in the elastic constants of other solids. 

First, we compare classical to quantum mechanical cal- 
culations in Fig. 0, where results in the NpT and the 
NVT ensemble are shown. In the simulations of Argon, 
a Lennard Jones potential V — 4e[(o'/r)^^ — [a/r)^] was 
used with parameters e = 1.67x 10^^^ J and a — 3.405 A. 
All simulations were based on system sizes TV = 500 and 
statistics of 5 x 10^ time steps. The Trotter number in the 
quantum runs has been chosen such that PT — 120. In- 
creasing P any further does not change the results within 
the statistical error bars. 
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FIG. 2. Elastic constant Cn of solid Argon at ambient 
pressure as a function of temperatre T. Classical and quan- 
tum simulations were carried out at constant volume V and 
constant external stress. 



In Fig. g, it can be seen that there is a considerable dis- 
crepency between classical and quantum mechanical elas- 
tic constants. The quantum Cn levels at about T = 15 K 
and one may extrapolate Cn « 3.5 GPa for the quantum 
mechanical ground state, while classically, Cn « 4.2 GPa 
at zero temperature. The relative effects in C12 and C44 
are similar. Note that this effect of a C(20%) reduction 
in Cii is considerably larger than the increase in the lat- 
tice constant of about 1% or the decrease of 10% of the 
heat of formation. 

The bulk of the reduction in Cn does not stem from 
the expressions introduced in Sec. [I A. There are two 



important contributions, namely the Born term and the 
term related to the fluctuation of the Lennard Jones po- 
tential. Details are given in Table |l|. The classical sys- 
tem at T = K has only one non-zero contribution, 
namely C\^^° — 4.2 GPa. As expected, differences be- 
tween classical and quantum results become negligible 
as the temperature approaches (or surpasses) the Debye 
temperature, which in the case of Argon is To — 93 K. 



B. Helium 



While Argon shows significant differences between clas- 
sical and quantum mechanical elastic constants, the cor- 
rections due to the kinetic energies are farily small. He- 
lium, in particular "^He, is much more "quantum" than 
Argon and therefore a more challenging test case than 
Argon. Xhe solid forms of ^He are only stable under 
pressure.ta Hence, we have to distinguish between elas- 
tic constants and Birch coefficients. There are three dif- 
ferent, stable '^He lattice structures: the bcc phase is 
mostly stable in the interval < T < 2 K at pressures 
3 MPa < p < 10 MPa, hep is the stable low-temperature 
phase at pressures p > 10 MPa, and there is a small tem- 
perature regime at pressures p > 150 MPa where a stable 
fee phase is found in between the hep and the fluid phase. 

The sipiulations for helium are all based on the Aziz 
potentiaO, which is known to be fairly reliable up to 
moderate pressures. Exchange effects are neglected. In 
the temperature regime accessible to our simulation they 
are generally believed to be unimportant. The transi- 
tion to a long-range ferromagnetic and antiferromagnetic 
transitions in hep and bcc '^He, respectively, only take 
place at temperatures in the mK regime. 



1. fee ^He 

In this subsection, we will focus on one particular rep- 
resentative point in the stable fee phase, because it is 
computationally very demanding to calculate elastic con- 
stants. In order to be able to work with relatively small 
Trotter numbers, it is necessary to go to large tempera- 
tures and low pressures. Yet, we want to avoid to be too 



close to a phase transition, e.g., the fee fluid phase tran- 
sition. The combination of T = 18 K and p = 200 MPa 
seems to be an appropriate choice: The quantum hmit is 
basicafly reached with Trotter number P = 32 and the 
transition to the fluid phase takes place at sufRciently 
higher temperature, namely at T = 22 K. Although the 
thermodynamic stability field of the hep phase is located 
at temperatures T < 18 K for the applied pressure, one 
may certainly expect metastability of the fee phase in the 
time window accesible to our simulations. 

The different Cy alo ng w ith their individual contribu- 



tions are listed in Table III . It is interesting to note that 
the Cauchy relation C12 — C44 for cubic crystals with 
central potentials basically also hold for the quantum 
solid, e.g., C44/C12 « 0.95, while the Birch coefficients 
don't: B44/B12 « 0.50. See Sec. |l^ for the definition 



and the relevance of Birch coefficients. It is instructive to 
represent the individual contributions to Cn graphically, 
which is done in Fig. H. It is noticable that the sum of 
the corrections to Cn which are related to the kinetic 
energy is relatively small, while the individual contribu- 
tions are fairly large. Yet, the solid is far away from being 
classical. The kinetic energy of ^He is about 88.9 fc^K, 
which is considerably larger than the thermal classical 
energy of 1.5 ks x 18 K = 27 ksK, thus T/Tr,ehye < 0.3. 
"Classical" helium would have Cn = 1.63 Gpa at the 
same external temperature and pressure. From Fig. 
we can see that it is possible to approximate the elastic 
constants fairly reasonably if only the Born contribution 
and the contribution due to the fiuctuating real potential 
are included into the calculation. 
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FIG. 3. Individual contributions to Cn of fee ^He at pres- 
sure p — 200 MPa and temperature T = 18 K. a: Born, b: 
flue, c: kin, d: Born-q, e: fluc-q, f: cross; corr. summarizes c, 
d, e, and f. Trotter Number P = 32 and number of particles 
N = 500. 



2. bcc ■^He 

The bcc phase can be considered to be the most in- 
teresting solid helium phase, because it is classically un- 
stable, fee and hep phase can both be considered to be 
classically stable at low temperatures, because their free 
energy differences are small. An interesting phenomena 



in the (quantum mechanical) bcc phase is the radial dis- 
tribution function g{r): The peaks in g(r) can not easily 
be related to nearest neighbors, next neighbors, etc., but 
the contributions of different neighboring shells "group" 
together. This is shown in Fig. 0. The final g{r) can 
be interpreted as a sum of broadened individual lines, 
which are represented in Fig. ^ as well. The overlap of 
such broadened lines is accompanied by a strong diffusion 
of individual atoms. We did not investigate in depth this 
diffusion process. 
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FIG. 4. Radial distribution function g{r) times r^ as 
a function of distance r for bcc "^He at T = 2 K and 
P = 10 MPa. The broken lines indicate the distance of near- 
est neighbors, next neighbors, etc. The width of the broken 
lines is proportional to the number of atoms in each shell. 

In the case of '^He, some experimental data is avail- 
able for the elastic constants or Birch coefficients. In 
Fig. ^, we compare the Bij calculated in our study 
with the available experimental data and some theoret- 
ical predicitions. Simulations at three different combi- 
nations of pressures and temperatures were performed: 
pi = 10 MPa with Ti = 2 K, p2 = 7 MPa with 
T2 = 1.5 K, and pa = 4 MPa with T3 = 1 K. In all 
cases. Trotter numbers of P = 256 were found to reffect 
the quantum limit sufficiently well, and the particle num- 
ber was TV = 432. The elastic constants can be assumed 
to be mainly temperature independent, as the temper- 
atures are far below the Debye temperatures, e.g., the 
ratios qi = (rkin)/l-5fcBr which are close to unity at 
the Debye temperature turned out to be: qi = 10.27, 
92 == 12.22, and q^ = 15.4. 

For bcc ■^He, the (relative) corrections in Cij which are 
related to the kinetic energy are much larger than in fee 
■^He. This can be seen in Table IV: the corrections make 
up nearly 80% of the total Cn for a pressure p = 4 MPa 
and about 30% for a pressure of p = lOPa (not explicitly 
shown in Tables). We want to note that the elastic con- 
stants obtained in the NVT ensemble and in the NpT 
ensemble agreed within the statistical error bars. 



3. hep ^He 

hep ^He has more independent elastic constants than 
bcc and fee ^He, respectively. The trend in all elastic 
constants, however, is yet again the same as in bcc or fee 
helium: the larger the pressure the more dominant the 
"classical" contribution to the elastic constants. Note 
that on the other hand, the absolute corrections due to 
the kinetic energy increase with increasing pressure. De- 
tails of the calculations are given in Tables M-^E Note 
that the classical elastic constants would be nearly three 
times as large as the quantum mechanical elastic con- 
stants at a pressure p = 90 MPa and even more than 
seven times as large in the case of p = 15 MPa. Due to 
the good agreement with the experimentally measured 
bulk modulus, see Fig. H, one may expect our quantum 
mechanical data to be fairly accurate. To our knowledge, 
the values presented in Tables ^- VI are predictions. No 
theoretical or experimental data is known to us. 
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FIG. 5. Birch coefficients for bcc 
and bulk modulus for hep 
tion of the molar volume, 
work. The broken line represents experimental data from 
Ref. [20] .The solid lines represent theoretical predictions 
from Ref. [21]. Open symbols refer to experimental data from 
Ref. [22]. 



C. Kinetic energy of solid He 

The calculation and measurment of kinetic energies 
(Tkin) in condensed helium phases has attracted a lot of 
recent attention. Condensed helium is highly quantum 
mechanical and therefore provides an ideal test ground 
for the application of many-body quantum statistical the- 
ories. Recent experimental developments have made it 
possible to measure kinetic energies with great accuracy 
and thus provide important tests on the theories. 

The data for (Tkin) presented in this study comple- 
ment data by Draeger and Ceperley, which has been 
presented recentlyE3. In particular, data for the bcc 
phase and the low-pressure hep regime are presented. 



see Fig. |[ We would like to emphasize that in both 
studies relatively high temperatures were taken in the 
high-pressure regime with small average atomic volume 
oi V < 20 A^, e.g., the ratio of kinetic energy to ther- 
mal energy, q — 2(rkin)/3fcBT is only slightly larger than 
three. Of course, this is still well below the Debye tem- 
perature, but some thermal activation will be present. 
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FIG. 6. Thermal kinetic energy (Tkin) as a function of the 
average volume per atom. Open diamonds are values from 
PIMC simulations by Draeger and Ceperley. The straight 
line is a power law fit to our data. 



Our data can be very well fit with an expression of the 
type 



(Tki 



AV 



(24) 



where A and a are fitparameters. A and a turn out to be 
a Ki 1.78 and A = 16100 if V is expressed in A^ per atom 
and kinetic energies in units of k^K. Within the regime 
considered here, our data is reflected within 1.5% accu- 
racy. In the bcc phase, this fit underestimated (Tkin) by 
about 1.5%, in the hep phase, (Tkin) is slightly overesti- 
mated. The value of (Tkin) near the quantum mechanical 
ground state seems to be mainly a function of the mo- 
lar or atomic volume only - relatively insensitive to the 
actual crystalline phase. 

We want to emphasize that the data points in Fig. g 
which turn out to be larger than the values suggested by 
the fit, all have relatively small values of the quantum 
parameter q. Thus, going to even smaller values of q, 
which is computationally expensive, might even increase 
the quality of the fit. Omitting all data with g < 5, an 
exponent of a = 1.75 is found. It is quite surprising that 
the birch coefficients depend exponentially on the molar 
volume (Fig. |^), the quantum mechanical ground state 
kinetic energy, however, only changes algebraicly with 
V. 



IV. CONCLUSIONS 

Expressions for the elastic constants of quantum solids 
have been presented in terms of a path integral represen- 
tation in the TVFT-ensemble. These expressions have 
then been applied and proven useful in path integral 
molecular dynamics simulations of solid Argon and hep, 
fee, and bcc helium III. In the NpT ensemble, the classi- 
cal formula can be used which relates strain fluctuations 
with elastic constants or - in the case of non-zero external 
pressure - with Birch coefficients. With the exception of 
hep and bcc '^He, the NVT expressions were dominated 
by the terms which one knows from classical simulations: 
the Born term and the potential fluctuation term. In 
the case of bcc '^He, terms related to the kinetic energy 
dominate the elastic constants. 

The quantum mechanical motion of the particles shows 
a stronger effect on the elastic constants than one might 
expect. E.g., in the case of Argon at zero external pres- 
sure, the quantum mechanical Cn is reduced by about 
20%, while cohesion energy and lattice constants are 
only decreased by 10% and increased by 1%, respectively. 
Hence, in order to have truely accurate estimates for elas- 
tic constants from computer simulations, quantum effects 
need to be taken into consideration. 

While deriving the expressions for elastic constants, an 
improved primitive estimator for the kinetic energy has 
been proposed. The statistical uncertainties of this es- 
timator do not increase with increasing Trotter number. 
Due to short correlation times, the improved primitive 
estimator results in statistical error bars smaller but in 
the order of the virial estimator for highly quantum sys- 
tems such as solid "^He. However, unlike the virial estima- 
tor, the improved primitive estimator can only be used 
in path integral molecular dynamics, but not in Monte 
Carlo simulations. 
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TABLE I. Individual components of dj for the ideal gas in units of NIcbT/V in the path integral representation. 





Born 


flue 


kin 


Born-q 


fluc-q 


cross 


corr 


sum 


quantum 


3.72639 
0.00005 


-0.30566 
0.00156 


0.04337 
0.00001 


0.03278 
0.00001 


-0.05530 
0.00053 


0.05469 
0.00150 


0.07555 
0.00137 


3.49627 
0.00195 


class. 


4.03116 
0.00004 


-0.15535 
0.00107 


0.00278 
0.00001 











0.00278 
0.00001 


3.87860 
0.00107 



TABLE II. Individual contributions to Cii for Argon at T = 7.5 K for the quantum and the classical calculation. Lower 
rows give statistical incertaintics based on 500.000 molecular dynamics steps. The individual contributions are introduced 
in Eq. (11). The colon named "corr" summarizes all corrections involving contributions from the kinetic energy: kin, Born-q, 
fluc-q, and cross. 





Born 


flue 


kin 


Born-q 


fluc-q 


cross 


corr 


sum 


Cii 

C44 


1.63 
0.77 
0.77 


-0.59 
-0.24 
-0.23 


0.43 


0.21 


0.38 


0.19 


-0.74 
0.01 
-0.36 


0.07 
0.12 
0.05 


0.14 
0.13 
0.09 


1.18 
0.66 
0.63 



TABLE III. Individual contributions to dj for fee He at T = 18 K and pressure p = 200 MPa. The colon named "corr" 
summarizes same terms as in previos table. Statistical error bars of final dj axe smaller than 5%. 





Born 


flue 


kin 


Born-q 


fluc-q 


cross 


corr 


sum 


Cii 


87.00 


-79.15 


88.59 


83.27 


-152.96 


8.26 


27.16 


35.00 


C12 


35.19 


-29.47 








4.84 


16.44 


21.28 


27.00 


C44 


35.19 


-27.84 


44.28 


41.62 


-78.53 


0.61 


7.98 


15.33 



TABLE IV. Individual contributions to dj for bee ^He at T = 2 K and pressure p = 4 MPa. The colon named "corr" 
summarizes same terms as in previos table. Statistical error bars of final dj are smaller than 10%. 



p 


T/K 


Cn 


C12 


Cl3 


C33 


C44 


Ce,6 


(l/>/cm^^ 


0.09 
0.05 
0.015 


10.0 
5.0 

2.5 


0.701 
0.405 
0.118 


0.303 
0.172 
0.052 


0.193 
0.111 

0.034 


0.822 
0.492 
0.144 


0.206 
0.123 
0.039 


0.193 
0.108 
0.032 


13.22 
14.85 
18.77 



TABLE V. Independent elastic constants of hep-'^He as obtained in the NpT ensemble. The thermal expectation value of 
the molar volume {V) is inserted as well. Pressure and elastic constants are expressed in GPa. Statistical error bars of Cij are 
smaller than 5%. 



F/cm^ 


T/K 


Cn 


C12 


Cl3 


C33 


C44 


Can 


(P> 


13.22 
14.85 
18.77 


10.0 
5.0 
2.5 


0.700 
0.407 
0.135 


0.234 
0.183 
0.043 


0.208 
0.125 
0.039 


0.810 
0.447 
0.156 


0.169 
0.138 
0.036 


0.192 
0.120 
0.025 


0.090 
0.050 
0.015 



TABLE VI. Independent elastic constants of hep-' He as obtained in the NVT ensemble. The thermal expectation value for 
the pressure (p) is inserted as well, p and Cij are expressed in GPa. 



